/**
* @file arc.cpp
* August 1, 2008
* 
* @brief Implementation of the ARC, a specialized curve object.
*
* This class derives the interface "Curve", implementing all the
* specific features of a ARC geometric object.
* 
* @author Antonio Miranda (amiranda@tecgraf.puc-rio.br )
* @author Andre Pereira (brabo@tecgraf.puc-rio.br)
* 
* @see Curve
* @see CurveFactory
*/

// includes
//----------
#include "arc.h"
#include <math.h>
#include <string.h>
#include <stdio.h>

// -------------------------------------------------------------------------
// Local variables and symbols:
//
#ifndef PI
#define PI (4.0*atan(1.0))
#endif

#define ARC_DISC PI*0.03

//
/////////////////////////////////////////////////////////////////////
Arc::Arc()
: m_numAddPoint( 0 )
{
  m_type = ARC;
}

//
/////////////////////////////////////////////////////////////////////
Arc::~Arc()
{
}

//////////////////////////////////////////////////////////////////////////
void Arc::getName( char *name )
{
  strcpy (name, "arc");
}

//
/////////////////////////////////////////////////////////////////////
vector<Point3D>
Arc::getCurvePoints()
{
  vector<Point3D> points(3);
  points[0] = m_center;
  points[1] = m_p1;
  points[2] = m_p2;
//  points[3] = m_normal;
  return points;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::getSubdvision( int nd,
                    double ratio,
                    vector<Point3D>& subdp )
{
  int  i;
  double x, y;
  Point3D vec;

  // check
  if (m_numAddPoint < 3)
  {
    subdp.resize (0);
    return;
  }

  // set dimension of vector
  subdp.resize (nd+1);

  // parametric subdvision
  vector<double> parSubd;
  getParamSubdivion (nd, ratio, parSubd);

  // begin point
  subdp[0] = m_p1;

  // intermediary points
  for (i = 1; i < nd; ++i)
  {
    // points on the plane
    x = (m_ratio * cos (m_angle*parSubd[i]));
    y = (m_ratio * sin (m_angle*parSubd[i]));

    // vector from center to desired point
    vec.setX( m_u.getX()*x + m_v.getX()*y);
    vec.setY( m_u.getY()*x + m_v.getY()*y);
    vec.setZ( m_u.getZ()*x + m_v.getZ()*y);

    // desired point in correct position
    subdp[i] = m_center + vec;
  }

  // end point
  subdp[nd] = m_p2;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::getTangent( int nd,
                 double ratio,
                 vector<Point3D>& tang )
{
  int  i;
  double x, y;
  Point3D vec;

  // check
  if (m_numAddPoint < 3)
  {
    tang.resize (0);
    return;
  }

  // set dimension of vector
  tang.resize (nd+1);

  // parametric subdvision
  vector <double> parSubd;
  getParamSubdivion (nd, ratio, parSubd);

  // intermediary points
  for (i = 1; i < nd+1; ++i)
  {
    // tangent on the plane
    x = -1.0 * (m_ratio * sin (m_angle*parSubd[i]));
    y = +1.0 * (m_ratio * cos (m_angle*parSubd[i]));

    // desired tangent on Cartesian space
    vec.setX( m_u.getX()*x + m_v.getX()*y );
    vec.setY( m_u.getY()*x + m_v.getY()*y );
    vec.setZ( m_u.getZ()*x + m_v.getZ()*y );

    // normalize tangent vector
    vec.normalize();
    tang[i] = vec;
  }
}

//
/////////////////////////////////////////////////////////////////////
double
Arc::getLenght()
{
  if (m_numAddPoint < 3)
    return 0.0;

  return m_angle * m_ratio;
}

// WHY DO WE NEED TOLERANCE HERE????
/////////////////////////////////////////////////////////////////////
int 
Arc::clstPoint( Point3D _pt,
                double& _t, 
                Point3D& _pt_on,
                double _tol )
{
  // project point to the plane
  Point3D onplane;
  pointToPlane( _pt, onplane );

  // performance checks
  if( m_center.distance2( onplane ) <= 0 )
    return 0;

  // unit vector from center to onplane point
  Point3D tmp = onplane - m_center;
  Point3D unit;
  tmp.normalize( unit );

  // set _pt_on in the correct position on plane and considering the ratio 
  _pt_on = m_center + unit * m_ratio;

  // compute angle
  double x = m_u.dotProd( _pt_on - m_center );
  double y = m_v.dotProd( _pt_on - m_center );
  double angle = myAtan2( y, x );

  // checks if the point is on the arc segment
  if( angle > m_angle )  // if it is not on the segment
  {
    // then, the closest point is one of the limiting points
    double distToP1 = onplane.distance( m_p1 );
    double distToP2 = onplane.distance( m_p2 );
    if( distToP1 > distToP2 )
    {
      _pt_on = m_p2;
      _t = 1.0;
    }
    else
    {
      _pt_on = m_p1;
      _t = 0.0;
    }
  }
  else
  {
    _t = angle / m_angle;
  }
  return 1;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::get3dPoint( double t,
                 Point3D& pt,
                 Point3D& tg )
{
  double angle = t * m_angle;
  double x = (m_ratio * cos (angle));
  double y = (m_ratio * sin (angle));

  // check
  if (m_numAddPoint < 3)
    return;
  
  // the point is p1 
  if (t == 0)
    pt = m_p1;

  // the point is p2 
  if (t == 1.0)
    pt = m_p2;

  // intermediary position  
  if (t> 0.0 && t < 1.0)
  {
    pt.setX( m_u.getX()*x + m_v.getX()*y );
    pt.setY( m_u.getY()*x + m_v.getY()*y );
    pt.setZ( m_u.getZ()*x + m_v.getZ()*y );
    pt = m_center + pt;
  }

  // local tangent vector
  double vx = -y;
  double vy =  x;
  tg.setX( m_u.getX()*vx + m_v.getX()*vy );
  tg.setY( m_u.getY()*vx + m_v.getY()*vy );
  tg.setZ( m_u.getZ()*vx + m_v.getZ()*vy );
  tg.normalize();
  tg *= this->getLenght();
}


//
/////////////////////////////////////////////////////////////////////
Point3D
Arc::get3dPoint( double _t )
{
  double angle = _t * m_angle;
  double x = (m_ratio * cos (angle));
  double y = (m_ratio * sin (angle));

  Point3D pt;

  // check
  if (m_numAddPoint < 3)
    return pt;

  // the point is p1 
  if (_t == 0)
    pt = m_p1;

  // the point is p2 
  if (_t == 1.0)
    pt = m_p2;

  // intermediary position  
  if (_t > 0.0 && _t < 1.0)
  {
    pt.setX( m_u.getX()*x + m_v.getX()*y );
    pt.setY( m_u.getY()*x + m_v.getY()*y );
    pt.setZ( m_u.getZ()*x + m_v.getZ()*y );
    pt = m_center + pt;
  }

  return pt;
}


//
/////////////////////////////////////////////////////////////////////
void 
Arc::draw()
{
  // check
  if (m_pDrawObj == NULL)
    return;

  // check
  if (m_numAddPoint < 3)
  {
    //std::cout << "_____np="<<m_numAddPoint << std::endl;
    return;
  }

  // draw arc
  drawArc(m_p2);
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::setBox()
{
  int i;
  m_box.clear();

  // discrete the arc
  double arc_disc = ARC_DISC;
  int    npts = (int) (m_angle / arc_disc);
  vector <Point3D> discPts;
  getSubdvision (npts, 1.0, discPts);

  // update to bound box 
  for (i = 0; i < npts+1; ++i)
     m_box.update( discPts[i].getX(), discPts[i].getY(), discPts[i].getZ() );
}




//
/////////////////////////////////////////////////////////////////////
int 
Arc::beginAddPoint( const vector<Point3D>& addp )
{
  if( addp.size() == 0 )
    return 0;

  addp[0].normalize( m_normal );

  m_numAddPoint = 0;
  return 1;
}

//
/////////////////////////////////////////////////////////////////////
int
Arc::addPoint( Point3D& pt )
{
  if( m_numAddPoint == 0 )
  {
    // checks if the normal was set, or if it is different of (0,0,0)
    if( m_normal.len() == 0.0)
      return 0;
    // adds the center point
    m_center = pt;
    // compute "d" plane parameter
    m_d = -1.0 * m_normal.dotProd( m_center );
  }
  else if( m_numAddPoint == 1 )
  {
    // set p1, but project point to the plane
    Point3D onplane;
    pointToPlane (pt, onplane);

    // performance checks
    m_ratio = m_center.distance2(onplane);
    if (m_ratio <= 0)
      return 0;
    else
      m_ratio = sqrt (m_ratio);

    m_p1 = onplane;

    // create local axis
    this->createLocalAxis();
    m_firstDraw = true;
  }
  else if( m_numAddPoint == 2 )
  {
    // project point to the plane
    Point3D onplane;
    pointToPlane( pt, onplane );

    // performance checks
    if (m_center.distance2(onplane) <= 0)
      return 0;

    // unit vector from center to onplane point
    Point3D tmp = onplane - m_center;
    Point3D unit;
    tmp.normalize(unit);

    // set p2 in the correct position on plane and considering the ratio 
    m_p2 = m_center + unit * m_ratio;

    // compute angle
    double x = m_u.dotProd( m_p2 - m_center );
    double y = m_v.dotProd( m_p2 - m_center );
    m_angle = myAtan2(y, x);

    if (m_angle == 0.0)
      m_angle = 2.0*PI;
    
    m_firstDraw = true;
  }
  else
    return 0;

  ++m_numAddPoint;
  
  return 1;
}

// TODO !!!
/////////////////////////////////////////////////////////////////////
int
Arc::addPoint( Point3D& pt,
               double tol )
{
  if( m_numAddPoint == 0 )
  {
    // checks if the normal was set, or if it is different of (0,0,0)
    if( m_normal.len() == 0.0)
      return 0;

    // adds the center point
    m_center = pt;
  
	// compute "d" plane parameter
    m_d = -1.0 * m_normal.dotProd( m_center );
  }
  else if( m_numAddPoint == 1 )
  {
    // set p1, but project point to the plane
    Point3D onplane;
    pointToPlane (pt, onplane);

    // performance checks
    m_ratio = m_center.distance2(onplane);

    if (m_ratio <= 0)
      return 0;
    else
      m_ratio = sqrt(m_ratio);

    m_p1 = pt;

    // create local axis
    Point3D tmp = m_p1 - m_center;
    tmp.normalize( m_u );  // u axes

    m_v = m_normal.crossProd( m_u );
    
	m_firstDraw = true;
  }
  else if( m_numAddPoint == 2 )
  {
    // set p2 in the correct position on plane and considering the ratio 
    m_p2 = pt;

    // compute angle
    double x = m_u.dotProd( m_p2 - m_center );
    double y = m_v.dotProd( m_p2 - m_center );

    m_angle = myAtan2(y,x);

    if( m_angle == 0.0 )
        m_angle = 2.0*PI;
    
    m_firstDraw = true;
  }
  else
    return 0;

  ++m_numAddPoint;
  
  return 1;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::addTempPoint( Point3D& pt )
{
  // checks
  if (m_numAddPoint < 1)
    return;
  if (m_pDrawObj == NULL)
    return;

  static int last = m_numAddPoint;

  int j;
  for (j = 0; j < 2; j++)
  {
    m_pDrawObj->setXor(1);

    if (m_firstDraw)
    {
      // draw with the center point
      if (m_numAddPoint == 1)
        drawFirst (pt);
      else if (m_numAddPoint == 2)
        drawSecond (pt);
      m_lastPt = pt;
      m_firstDraw = false;
      m_pDrawObj->end();
      m_pDrawObj->setXor(0);
      break;
    }
    else
    {
      if (j == 0)
      {
        // draw with the center point
        if (m_numAddPoint == 1)
          drawFirst (m_lastPt);
        else if (m_numAddPoint == 2)
          drawSecond (m_lastPt);
      }
      else
      {
        // draw with the center point
        if (m_numAddPoint == 1)
          drawFirst (pt);
        else if (m_numAddPoint == 2)
          drawSecond (pt);
        m_lastPt = pt;
      }
    }
    m_pDrawObj->setXor(0);
  }
}


//
/////////////////////////////////////////////////////////////////////
void
Arc::addTempOrthoPoint( Point3D& pt )
{
  // checks
  if (m_numAddPoint < 1)
    return;
  if (m_pDrawObj == NULL)
    return;

  static int last = m_numAddPoint;

  int j;
  for (j = 0; j < 2; j++)
  {
    m_pDrawObj->setXor(1);

    if (m_firstDraw)
    {
      // draw with the center point
      if (m_numAddPoint == 1)
        drawFirst (pt);
      else if (m_numAddPoint == 2)
        drawSecond (pt);
      m_lastPt = pt;
      m_firstDraw = false;
      m_pDrawObj->end();
      m_pDrawObj->setXor(0);
      break;
    }
    else
    {
      if (j == 0)
      {
        // draw with the center point
        if (m_numAddPoint == 1)
          drawFirst (m_lastPt);
        else if (m_numAddPoint == 2)
          drawSecond (m_lastPt);
      }
      else
      {
        Point3D pt_diff = pt - m_center;

        double abs_x = fabs(pt_diff.getX());
        double abs_y = fabs(pt_diff.getY());
        double abs_z = fabs(pt_diff.getZ());

        if ( (abs_x > abs_y) && (abs_x > abs_z) )
        {
          pt_diff.setY(0);
          pt_diff.setZ(0);
        }
        else if ( abs_y > abs_z )
        {
          pt_diff.setX(0);
          pt_diff.setZ(0);
        }
        else
        {
          pt_diff.setX(0);
          pt_diff.setY(0);
        }

        pt = m_center + pt_diff;

        // draw with the center point
        if (m_numAddPoint == 1)
          drawFirst (pt);
        else if (m_numAddPoint == 2)
          drawSecond (pt);
        m_lastPt = pt;
      }
    }
    m_pDrawObj->setXor(0);
  }
}

//
/////////////////////////////////////////////////////////////////////
Point3D
Arc::getLastTempPoint()
{ 
  return m_lastPt;
}

// end add point and the curve is finally created
/////////////////////////////////////////////////////////////////////
int
Arc::endAddPoint()
{
  if( !possible() )
    return 0;
  setBox();
  return 1;
}

/////////////////////////////////////////////////////////////////////
void
Arc::getEndPoints( Point3D pt[2] )
{
  pt[0] = m_p1;
  pt[1] = m_p2;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::delPoint( Point3D& pt )
{
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::delPoint( Point3D& pt,
               double tol )
{
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::copy( Curve** cp )
{
  *cp = new Arc;
  Arc *subobj = (Arc *) *cp;
  subobj->m_normal = this->m_normal;
  subobj->m_center = this->m_center;
  subobj->m_p1 = this->m_p1; 
  subobj->m_p2 = this->m_p2;
  subobj->m_ratio = this->m_ratio;
  subobj->m_angle = this->m_angle;
  subobj->m_normal = this->m_normal;
  subobj->m_d = this->m_d;
  subobj->m_u = this->m_u;
  subobj->m_v = this->m_v;
  subobj->m_numAddPoint = this->m_numAddPoint;
  subobj->m_type = this->m_type;
  subobj->m_subdv = this->m_subdv;
  subobj->m_ratio = this->m_ratio;
  subobj->m_attr = this->m_attr;
  subobj->m_box = this->m_box;
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::transform( double m[4][4] )
{
  Point3D c  = m_center;
  Point3D pd = c;
  pd += m_normal;

  Point3D a = pd;
  double w = m[3][0]*a.getX()+m[3][1]*a.getY()+m[3][2]*a.getZ()+m[3][3];
  pd.setX ((m[0][0]*a.getX() + m[0][1]*a.getY() + m[0][2]*a.getZ() + m[0][3]) / w);
  pd.setY ((m[1][0]*a.getX() + m[1][1]*a.getY() + m[1][2]*a.getZ() + m[1][3]) / w);
  pd.setZ ((m[2][0]*a.getX() + m[2][1]*a.getY() + m[2][2]*a.getZ() + m[2][3]) / w);
  a = c;

  w = m[3][0]*a.getX()+m[3][1]*a.getY()+m[3][2]*a.getZ()+m[3][3];
  c.setX ((m[0][0]*a.getX() + m[0][1]*a.getY() + m[0][2]*a.getZ() + m[0][3]) / w);
  c.setY ((m[1][0]*a.getX() + m[1][1]*a.getY() + m[1][2]*a.getZ() + m[1][3]) / w);
  c.setZ ((m[2][0]*a.getX() + m[2][1]*a.getY() + m[2][2]*a.getZ() + m[2][3]) / w);
  m_normal = pd - c;
  m_normal.normalize();

  a = m_center;
  w = m[3][0]*a.getX()+m[3][1]*a.getY()+m[3][2]*a.getZ()+m[3][3];
  m_center.setX ((m[0][0]*a.getX() + m[0][1]*a.getY() + m[0][2]*a.getZ() + m[0][3]) / w);
  m_center.setY ((m[1][0]*a.getX() + m[1][1]*a.getY() + m[1][2]*a.getZ() + m[1][3]) / w);
  m_center.setZ ((m[2][0]*a.getX() + m[2][1]*a.getY() + m[2][2]*a.getZ() + m[2][3]) / w);

  a = m_p1;
  w = m[3][0]*a.getX()+m[3][1]*a.getY()+m[3][2]*a.getZ()+m[3][3];
  m_p1.setX ((m[0][0]*a.getX() + m[0][1]*a.getY() + m[0][2]*a.getZ() + m[0][3]) / w);
  m_p1.setY ((m[1][0]*a.getX() + m[1][1]*a.getY() + m[1][2]*a.getZ() + m[1][3]) / w);
  m_p1.setZ ((m[2][0]*a.getX() + m[2][1]*a.getY() + m[2][2]*a.getZ() + m[2][3]) / w);

  a = m_p2;
  w = m[3][0]*a.getX()+m[3][1]*a.getY()+m[3][2]*a.getZ()+m[3][3];
  m_p2.setX ((m[0][0]*a.getX() + m[0][1]*a.getY() + m[0][2]*a.getZ() + m[0][3]) / w);
  m_p2.setY ((m[1][0]*a.getX() + m[1][1]*a.getY() + m[1][2]*a.getZ() + m[1][3]) / w);
  m_p2.setZ ((m[2][0]*a.getX() + m[2][1]*a.getY() + m[2][2]*a.getZ() + m[2][3]) / w);
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::invert()
{
  // invert normal
  m_normal *= -1.0;
  m_d *= -1.0;

  // swap points
  Point3D swap = m_p1;
  m_p1 = m_p2;
  m_p2 = swap;

  this->createLocalAxis();
}

//
/////////////////////////////////////////////////////////////////////
int
Arc::possible()
{
   //WARNING add this check: m_dir.normalize() 
  if( m_numAddPoint == 3 )
    return 1;

  return 0;
}

//
/////////////////////////////////////////////////////////////////////
void
Arc::print()
{
  std::cout << " - Arc: " << std::endl;
  std::cout << "   n: ( " << m_normal.getX() << "; " << m_normal.getY() << "; " << m_normal.getZ() << " )" << std::endl;
  std::cout << "   c: ( " << m_center.getX() << "; " << m_center.getY() << "; " << m_center.getZ() << " )" << std::endl;
  std::cout << "   1: ( " << m_p1.getX()     << "; " << m_p1.getY()     << "; " << m_p1.getZ()     << " )" << std::endl;
  std::cout << "   2: ( " << m_p2.getX()     << "; " << m_p2.getY()     << "; " << m_p2.getZ()     << " )" << std::endl;
  std::cout << "   u: ( " << m_u.getX()      << "; " << m_u.getY()      << "; " << m_u.getZ()      << " )" << std::endl;
  std::cout << "   v: ( " << m_v.getX()      << "; " << m_v.getY()      << "; " << m_v.getZ()      << " )" << std::endl;
  std::cout << "   r: " << m_ratio << std::endl;
  std::cout << "   a: " << m_angle << std::endl;
}




//
/////////////////////////////////////////////////////////////////////
int 
Arc::planeIntersect( Point3D& _n,
                    double _d,
                    vector<Curve*>& _curves,
                    vector<Point3D>& _pts)
{
  // using algorithms: "The intersection of two planes" and
  // "The intersection of a line and a sphere (or a circle)"
  // from http://local.wasp.uwa.edu.au/~pbourke/geometry/

  // return statements:
  // 0 : the segment lies strictly to one side or the other of the plane.
  // 1 : the segment intersects the plane only once (or the plane is tangent to the arc).
  // 2 : the segment intersects the plane twice.
  // 3 : the segment lies wholly within the plane.
  
  _curves.clear();
  _n.normalize();

  //-------------------------------------------------
  // first: finds the intersection of two planes
  //        N1 . p = d1    and    N2 . p = d2
  //-------------------------------------------------
  // defines the direction of the intersection line (N1 * N2)
  Point3D dirL = m_normal.crossProd( _n );
  // checks if the planes are parallel or coincident (N1 * N2 = 0)
  if( dirL.len() == 0.0 )  // if parallel
  {
    if( m_center.dotProd( _n ) == _d ) // if coincident
      return 3;
    else
      return 0;
  }
  // construct the LINE equation: p(t) = c1 N1 + c2 N2 + t (N1 * N2)
  // to finds a point, do a linear combination with the two normal vectors
  double N1dotN1 = m_normal.dotProd( m_normal );
  double N2dotN2 = _n.dotProd( _n );
  double N1dotN2 = m_normal.dotProd( _n );
  double det = N1dotN1 * N2dotN2 - N1dotN2 * N1dotN2;
  double c1 = ( m_d * N2dotN2 - _d  * N1dotN2 ) / det;
  double c2 = ( _d  * N1dotN1 - m_d * N1dotN2 ) / det;
  Point3D p0( c1 * m_normal + c2 * _n );

  //-----------------------------------------------------------------
  // second: finds the intersection of a line and a sphere (or circle)
  //         p(t)=p1+t(p2-p1)  and  (x-x3)^2+(y-y3)^2+(y-y3)^2=r^2
  //-----------------------------------------------------------------
  // combining the two equations, it results in:  a*t^2 + b*t + c = 0
  // therefore, one should solve the quadratic equation and study its behaviour
  double a = pow( dirL.len(), 2 );
  double b = 2 * dirL.dotProd( p0 - m_center );
  double c = m_center.dotProd( m_center ) + p0.dotProd( p0 ) - pow( m_ratio, 2 )
           - 2 * p0.dotProd( m_center );
  double delta = b * b - 4.0 * a * c;

  if( delta < 0.0 )  // the plane does not intersect the circle
    return 0;
  else if( delta == 0.0 )  // the plane is tangent to the circle <- THE PLACE TO ADD TOLERANCE
  {
    double t = -b/2.0/a;
    Point3D p1( p0 + t * dirL );
    // verifies, finally, if the point lies on the arc
    // compute angle
    double x = m_u.dotProd( p1 - m_center );
    double y = m_v.dotProd( p1 - m_center );
    double angle = myAtan2( y, x );
    if( angle > m_angle )  // the point is outside
      return 0;
    else
    {
      // creates two new curves (ARC) and adds the curves to the list
      Arc* curve1 = this->createArc( m_p1, p1 );
      Arc* curve2 = this->createArc( p1, m_p2 );
      _curves.push_back( curve1 );
      _curves.push_back( curve2 );
      return 1;
    }
  }
  else  // the plane intersect twice the circle
  {
    double t1 = (-b+sqrt(delta))/2/a;
    double t2 = (-b-sqrt(delta))/2/a;
    Point3D p1( p0 + t1 * dirL );
    Point3D p2( p0 + t2 * dirL );
    // verifies, finally, if the point lies on the arc, by computing the angle
    double x1 = m_u.dotProd( p1 - m_center );
    double y1 = m_v.dotProd( p1 - m_center );
    double angle1 = myAtan2( y1, x1 );
    double x2 = m_u.dotProd( p2 - m_center );
    double y2 = m_v.dotProd( p2 - m_center );
    double angle2 = myAtan2( y2, x2 );
    if( (angle1 > m_angle) && (angle2 > m_angle) )  // the point is outside
      return 0;
    else if( angle1 > m_angle )  // only the point "p1" is outside
    {
      // creates two new curves (ARC) and adds the curves to the list
      Arc* curve1 = this->createArc( m_p1, p2 );
      Arc* curve2 = this->createArc( p2, m_p2 );
      _curves.push_back( curve1 );
      _curves.push_back( curve2 );
      return 1;
    }
    else if( angle2 > m_angle )  // only the point "p2" is outside
    {
      // creates two new curves (ARC) and adds the curves to the list
      Arc* curve1 = this->createArc( m_p1, p1 );
      Arc* curve2 = this->createArc( p1, m_p2 );
      _curves.push_back( curve1 );
      _curves.push_back( curve2 );
      return 1;
    }
    else
    {
      // creates three new curves (ARC) and adds the curves to the list
      Arc* curve1 = NULL;
      Arc* curve2 = NULL;
      Arc* curve3 = NULL;
      if( angle1 > angle2 )
      {
        curve1 = this->createArc( m_p1, p2 );
        curve2 = this->createArc( p2, p1 );
        curve3 = this->createArc( p1, m_p2 );
      }
      else
      {
        curve1 = this->createArc( m_p1, p1 );
        curve2 = this->createArc( p1, p2 );
        curve3 = this->createArc( p2, m_p2 );
      }
      _curves.push_back( curve1 );
      _curves.push_back( curve2 );
      _curves.push_back( curve3 );
      return 2;
    }
  }
}

//
/////////////////////////////////////////////////////////////////////
int
Arc::planeIntersect( Point3D& _orig,
                    Point3D& _n,
                    vector<Curve*>& _curves,
                    vector<Point3D>& _pts)
{
  // finds the minimal distance from origin to the plane
  _n.normalize();
  double d = _orig.dotProd( _n );
  return planeIntersect( _n, d, _curves, _pts);
}

//
/////////////////////////////////////////////////////////////////////
int 
Arc::join( Curve* _givenCurve,
           Curve** _thirdCurve )
{
  return 1;
}

//
/////////////////////////////////////////////////////////////////////
Curve*
Arc::join( Curve* _otherCurve, double _tol )
{
  Point3D pi;
  Point3D ps;
  Point3D pe; 
  Arc* arc = NULL;

  Arc* c2 = (Arc *)_otherCurve;

  if( m_normal.distance( c2->m_normal ) > _tol )
    return 0;

  if( m_center.distance( c2->m_center ) > _tol )
    return 0;

  if( m_p1.isAtSamePosition(c2->m_p1, _tol) )
  {
    pi = m_p1;
    ps = m_p2;
    pe = c2->m_p2;
  }
  else if( m_p1.isAtSamePosition(c2->m_p2, _tol) )
  {
    pi = m_p1;
    ps = m_p2;
    pe = c2->m_p1;
  }
  else if( m_p2.isAtSamePosition(c2->m_p1, _tol) )
  {
    pi = m_p2;
    ps = m_p1;
    pe = c2->m_p2;
  }
  else if( m_p2.isAtSamePosition(c2->m_p2, _tol) )
  {
    pi = m_p2;
    ps = m_p1;
    pe = c2->m_p1;
  }
  else
  {
    std::cout << "ARC::join() -> merge of arc not possible!!" << std::endl;
    return NULL;
  }

  // closed curve. Nothing...
  if( ps == pe )
    return 0;

  // create local axis
  Point3D u = ps - m_center;
  u.normalize();  // u axes
  Point3D v = m_normal.crossProd( u );
  // compute angle
  double x1 = u.dotProd( pe - m_center );
  double y1 = v.dotProd( pe - m_center );
  double ang1 = myAtan2(y1, x1);
  double x2 = u.dotProd( pi - m_center );
  double y2 = v.dotProd( pi - m_center );
  double ang2 = myAtan2(y2, x2);

  if( ang2 < ang1 )
    arc = this->createArc( ps, pe );
  else
    arc = this->createArc( pe, ps );

  return arc;
}


// TODO: checks return statement, ONLY SPLIT if CURVE is completed?
/////////////////////////////////////////////////////////////////////
int 
Arc::split( vector<Curve*>& _curves )
{
  split( m_subdv, m_ratio, _curves );
  return 1;
}

// TODO: checks return statement, ONLY SPLIT if CURVE is completed?
/////////////////////////////////////////////////////////////////////
int 
Arc::split( int _subdvision,
            double _ratio,
            vector<Curve*>& _curves )
{
  // get the subdivision points for the curve
  vector<Point3D> subv;
  getSubdvision( _subdvision, _ratio, subv );
  
  // loop over the number of subdivision points to build the new curves
  int size = (int) subv.size();
  for( int i=1; i<size; i++ )
  {
    // creates a new curve (ARC) and adds the curve to the list
    Arc* curve = this->createArc( subv[i-1], subv[i] );
    _curves.push_back( curve );
  }
  return 1;
}

// TODO: checks return statement, ONLY SPLIT if CURVE is completed?
/////////////////////////////////////////////////////////////////////

Curve* 
Arc::split( double _t0, double _t1, double _tol )
{
  Curve* arc = NULL;

  // if the interval is from 0 to 1 then it creates a copy of the curve
  if( (_t0==0.0) && (_t1==1.0) )
  {
    this->copy( &arc );
  }
  else
  {
    // creates a new curve (ARC)
    arc = new Arc;

    // adds normal, center and two points
    arc->beginAddPoint( m_normal.getX(), m_normal.getY(), m_normal.getZ() );
    arc->addPoint( m_center );
    Point3D _pts =  get3dPoint(_t0);
    arc->addPoint( _pts );
    _pts =  get3dPoint(_t1);
    arc->addPoint( _pts );
    arc->endAddPoint();
  }

  return arc;
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::pointToPlane( const Point3D& pt,
                   Point3D& onplane )
{
  Point3D vec_pt = pt - m_center;

  // distance of point to the plane
  double dist = m_normal.dotProd(vec_pt);
  Point3D sub_vec = m_normal * dist;

  onplane = pt - sub_vec;
}

//
/////////////////////////////////////////////////////////////////////
double
Arc::myAtan2( double y,
              double x )
{
  double at = atan2(y,x);

  if (at < 0) at += 2*PI;

  return(at);
}

//
/////////////////////////////////////////////////////////////////////
void 
Arc::drawArc( Point3D &endPt )
{
  double x = m_u.dotProd(endPt - m_center); // point on local axis 
  double y = m_v.dotProd(endPt - m_center);
  double angle = myAtan2(y, x);         // current angle
  double arc_disc = ARC_DISC;

  if( m_p1.isAtSamePosition(m_p2,0.00001*ARC_DISC))
    angle = 2.0*PI;

  int npts = (int) (angle / arc_disc);  // discrete the arc

  // check of draw obj
  if (m_pDrawObj == NULL)
    return;

  if (npts == 0)
    npts = 1;
 
  double dang = angle / npts; // delta angles
  int    i;
  Point3D currpts;

  m_pDrawObj->begin();
  for (i = 0; i < npts; ++i)
  {
    // local coordinate position on local plane
    x = (m_ratio * cos (i*dang));
    y = (m_ratio * sin (i*dang));

    // vector from center to 3d position
    currpts.setX( m_u.getX()*x + m_v.getX()*y );
    currpts.setY( m_u.getY()*x + m_v.getY()*y );
    currpts.setZ( m_u.getZ()*x + m_v.getZ()*y );

    // add center to last vector
    currpts = m_center + currpts;
    
    // draw the segment
    m_pDrawObj->addVertex(currpts.getX(), currpts.getY(), currpts.getZ());
  }
  // end point
  m_pDrawObj->addVertex(endPt.getX(), endPt.getY(), endPt.getZ());
  
  m_pDrawObj->end();
}

//
//////////////////////////////////////////////////////////////////////////
void 
Arc::drawFirst( Point3D &pt )
{
  static char text[50];
  static double ratio;
  m_pDrawObj->begin();
  m_pDrawObj->addVertex(m_center.getX(), m_center.getY(), m_center.getZ());
  m_pDrawObj->addVertex(pt.getX(), pt.getY(), pt.getZ());
  m_pDrawObj->end();

  ratio = m_center.distance(pt);

  sprintf(text, "R=%f", ratio);

  Point3D c = (m_center + pt ) / 2.0;

  m_pDrawObj->text(c.getX(), c.getY(), c.getZ(), text);
}

//
//////////////////////////////////////////////////////////////////////////
void 
Arc::drawSecond( Point3D &pt )
{
  // project point to the plane
  Point3D onplane;
  pointToPlane (pt, onplane);

  // performance checks
  if (m_center.distance2(onplane) <= 0)
    return;

  // unit vector from center to onplane point
  Point3D tmp = onplane - m_center;
  Point3D unit;
  tmp.normalize(unit);

  // correct position on plane and considering the ratio 
  Point3D correctPt = m_center + unit * m_ratio;

  // draw text
  double x = m_u.dotProd( correctPt - m_center ); // point on local axis 
  double y = m_v.dotProd( correctPt - m_center );
  double angle = myAtan2(y, x)* 180.0 / PI;         // current angle
  static char text[50];
  sprintf(text, "ang=%.3f", angle);
  //printf ("ang=%f", angle);
  Point3D c = (m_center + correctPt) / 2.0;
  m_pDrawObj->text(c.getX(), c.getY(), c.getZ(), text);

  drawArc(correctPt);

  m_pDrawObj->begin();
  m_pDrawObj->addVertex(m_center.getX(), m_center.getY(), m_center.getZ());
  m_pDrawObj->addVertex(m_p1.getX(), m_p1.getY(), m_p1.getZ());
  m_pDrawObj->end();

  m_pDrawObj->begin();
  m_pDrawObj->addVertex(m_center.getX(), m_center.getY(), m_center.getZ());
  m_pDrawObj->addVertex(correctPt.getX(), correctPt.getY(), correctPt.getZ());
  m_pDrawObj->end();
}

//
//////////////////////////////////////////////////////////////////////////
Arc*
Arc::createArc( Point3D& _p1,
                Point3D& _p2 )
{
  // creates a new curve (ARC)
  Arc* arc = new Arc();
  // adds normal, center and two points
  arc->m_normal = this->m_normal;
  arc->addPoint( m_center );
  arc->addPoint( _p1 );
  arc->addPoint( _p2 );
  return arc;
}





//
/////////////////////////////////////////////////////////////////////
void
Arc::updateData( int dg, int nP, int nU, double* P, double* U )
{
  // setting the read data in the curve variables:
  if( 12 == nP ) // if it is not reading from Topology
  {
    m_center.setXYZ( P[0], P[1], P[2] );  
    m_p1.setXYZ( P[3], P[4], P[5] );     
    m_p2.setXYZ( P[6], P[7], P[8] );
    m_normal.setXYZ( P[9], P[10], P[11] );
  }
  else
  {
    m_p1.setXYZ( P[0], P[1], P[2] );     
    m_p2.setXYZ( P[4*(nP-1)], P[4*(nP-1)+1], P[4*(nP-1)+2] );
  }
}




vector<Point3D> 
Arc::getPolygonalPoints()
{ 
	vector<Point3D> temp;
	double x = m_u.dotProd(m_p2 - m_center); // point on local axis 
	double y = m_v.dotProd(m_p2 - m_center);
	double angle = myAtan2(y, x);         // current angle
	double arc_disc = ARC_DISC;
	if( m_p1.isAtSamePosition(m_p2,0.00001*ARC_DISC))
		angle = 2.0*PI;
	int npts = (int)(angle / arc_disc)+1;  // discrete the arc
	if (npts == 0)
		npts = 1;
	double dang = angle / (npts-1); // delta angles
	int    i;
	Point3D currpts;
	for (i = 0; i < npts; ++i)
	{
		// local coordinate position on local plane
		x = (m_ratio * cos (i*dang));
		y = (m_ratio * sin (i*dang));
		// vector from center to 3d position
		currpts.setX( m_u.getX()*x + m_v.getX()*y );
		currpts.setY( m_u.getY()*x + m_v.getY()*y );
		currpts.setZ( m_u.getZ()*x + m_v.getZ()*y );
		// add center to last vector
		currpts = m_center + currpts;
        temp.push_back(currpts);
	}
	return temp;
}

//
//////////////////////////////////////////////////////////////////////////
void Arc::createLocalAxis()
{
  Point3D tmp = m_p1 - m_center;
  tmp.normalize( m_u );  // u axes
  m_v = m_normal.crossProd( m_u );
}

